Spectroscopic Data Reduction#
Spectroscopic data reduction is the process of converting raw spectroscopic images from a spectrograph into a calibrated spectrum. This process involves several steps, including bias subtraction, flat fielding, wavelength calibration, and extraction of the spectrum. In this notebook, we will go through the process of reducing a spectroscopic image using Python.
The CCD pre-processing steps are as usual: bias/dark subtraction, flat fielding, and cosmic ray removal. The spectroscopic data reduction steps are: extraction of the spectrum, wavelength calibration and flux calibration. Here we will focus mostly on the spectroscopic data reduction steps. After the data reduction, the spectrum can be used for scientific analysis, such as measuring redshifts, line identifications, and measuring line fluxes. In this notebook, we will reduce a 2D spectrum of an afterglow of a gamma-ray burst (GRB). Since this tutorial is prepared for the Project No. 5 of our TAO class, I leave the scientific analysis of the spectrum to readers.
The data used in this notebook is from the Subaru Telescope’s Faint Object Camera and Spectrograph (FOCAS). The data is 2D images of a spectrum of an afterglow of a GRB, GRB 130606A, which is a significantly redshifted source at \(z\sim5.9\). The data was taken on June 7, 2013, using the VPH900 grism, which covers the wavelength range of 7500-10450 Angstroms (Subaru FOCAS Hompage).
Let’s start by organizing the data and displaying the 2D spectral images.
import warnings
from datetime import datetime
from pathlib import Path
import numpy as np
import preproc # custom module for data reduction
from astropy.io import ascii, fits
from astropy.modeling.fitting import LevMarLSQFitter
from astropy.modeling.models import Chebyshev2D, Gaussian1D
from astropy.nddata import CCDData
from astropy.stats import gaussian_fwhm_to_sigma, sigma_clip
from astropy.table import Column, Table
from astropy.visualization import ZScaleInterval
from ccdproc import combine, cosmicray_lacosmic, gain_correct
from IPython.display import Image
from matplotlib import gridspec
from matplotlib import pyplot as plt
from matplotlib import rc, rcParams
from mpl_toolkits.axes_grid1 import make_axes_locatable
from numpy.polynomial.chebyshev import chebfit, chebval
from scipy.interpolate import UnivariateSpline
from scipy.ndimage import median_filter
from scipy.optimize import curve_fit
from scipy.signal import peak_widths
from skimage.feature import peak_local_max
from astropy import units as u
# plt.rcParams['axes.linewidth'] = 2
plt.rcParams["font.size"] = 12
plt.rcParams["axes.labelsize"] = 14
WD = Path("./data/proj5")
RAWDIR = WD / "raw" / "P03hbahk0304145443GZ_data"
OUTDIR = WD / "out"
if not OUTDIR.exists():
OUTDIR.mkdir(parents=True)
warnings.filterwarnings("ignore", category=UserWarning, append=True)
def make_summary_table_focas(rawdir, suffix=".fits.gz"):
"""Make a summary table of the raw data in the directory. This function is
specifically designed for the data taken with the FOCAS instrument.
Args:
rawdir (pathlib.Path): The directory containing the raw data.
suffix (str, optional): The suffix of the raw data files. Defaults to
'.fits.gz'.
Returns:
stab (astropy.table.Table): The summary table of the raw data in the directory.
"""
# making a summary table
summary = []
for f in rawdir.glob("*" + suffix):
hdr = fits.getheader(f)
# getting the filter information from the header
filt1 = hdr["FILTER01"]
filt2 = hdr["FILTER02"]
filt3 = hdr["FILTER03"]
X = hdr["AIRMASS"] # airmass
disp = hdr["DISPERSR"] # disperser
# appending the data to the summary list
summary.append(
[
f.name,
hdr["DATE-OBS"],
hdr["OBJECT"],
hdr["RA"],
hdr["DEC"],
hdr["DATA-TYP"],
hdr["EXPTIME"],
X,
hdr["ALTITUDE"],
disp,
hdr["SLT-WID"],
hdr["SLT-PA"],
hdr["SLIT"],
filt1,
filt2,
filt3,
hdr["UT-STR"],
hdr["UT-END"],
]
)
# creating the summary table
stab = Table(
rows=summary,
names=[
"filename",
"date-obs",
"object",
"ra",
"dec",
"type",
"exptime",
"airmass",
"altitude",
"disperser",
"slit-width",
"slit-pa",
"slit",
"filter1",
"filter2",
"filter3",
"ut-str",
"ut-end",
],
dtype=[
"U50",
"U50",
"U50",
"U50",
"U50",
"U50",
"f8",
"f8",
"f8",
"U50",
"f8",
"f8",
"U50",
"U50",
"U50",
"U50",
"U50",
"U50",
],
)
return stab
summary_table = make_summary_table_focas(RAWDIR)
summary_table
| filename | date-obs | object | ra | dec | type | exptime | airmass | altitude | disperser | slit-width | slit-pa | slit | filter1 | filter2 | filter3 | ut-str | ut-end |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| str50 | str50 | str50 | str50 | str50 | str50 | float64 | float64 | float64 | str50 | float64 | float64 | str50 | str50 | str50 | str50 | str50 | str50 |
| FCSA00141794.fits.gz | 2013-06-07 | DOMEFLAT | 11:17:48.038 | +19:54:08.26 | DOMEFLAT | 3.0 | 1.0 | 89.9361 | SCFCGRHD90 | 0.5 | -90.3 | SCFCSLLC08 | NONE | SCFCFLSO58 | NONE | 04:37:11.906 | 04:37:15.139 |
| FCSA00141920.fits.gz | 2013-06-07 | GRB 130606A | 16:37:41.157 | +29:47:45.71 | OBJECT | 1200.0 | 1.136 | 61.6129 | SCFCGRHD90 | 0.5 | -90.3 | SCFCSLLC08 | NONE | SCFCFLSO58 | NONE | 07:59:00.676 | 08:19:01.117 |
| FCSA00141930.fits.gz | 2013-06-07 | GRB 130606A | 16:37:41.230 | +29:47:45.13 | OBJECT | 1200.0 | 1.038 | 74.4191 | SCFCGRHD90 | 0.5 | -90.3 | SCFCSLLC08 | NONE | SCFCFLSO58 | NONE | 09:03:28.980 | 09:23:29.405 |
| FCSA00141848.fits.gz | 2013-06-07 | BIAS | 11:39:50.841 | +19:54:11.18 | BIAS | 0.0 | 1.0 | 89.9361 | SCFCGRMR01 | 0.5 | -90.3 | SCFCSLLC20 | NONE | SCFCFLSO58 | NONE | 04:59:10.571 | 04:59:10.574 |
| FCSA00141916.fits.gz | 2013-06-07 | GRB 130606A | 16:37:41.020 | +29:47:44.75 | OBJECT | 1200.0 | 1.254 | 52.8249 | SCFCGRHD90 | 0.5 | -90.3 | SCFCSLLC08 | NONE | SCFCFLSO58 | NONE | 07:18:00.165 | 07:38:00.608 |
| FCSA00141734.fits.gz | 2013-06-06 | BIAS | 21:41:03.767 | +17:46:28.41 | BIAS | 0.0 | 1.001 | 86.8997 | SCFCGRHD90 | 0.5 | -0.3 | SCFCSLLC20 | NONE | SCFCFLSO58 | NONE | 14:53:03.807 | 14:53:03.810 |
| FCSA00141566.fits.gz | 2013-06-06 | DOMEFLAT | 11:23:03.374 | +19:54:09.25 | DOMEFLAT | 1.1 | 1.0 | 89.9625 | SCFCGRHD90 | 0.5 | -90.3 | SCFCSLLC20 | NONE | SCFCFLSO58 | NONE | 04:46:28.971 | 04:46:30.721 |
| FCSA00141576.fits.gz | 2013-06-06 | DOMEFLAT | 11:26:18.779 | +19:54:09.76 | DOMEFLAT | 3.0 | 1.0 | 89.9625 | SCFCGRHD90 | 0.5 | -90.3 | SCFCSLLC08 | NONE | SCFCFLSO58 | NONE | 04:49:43.639 | 04:49:46.993 |
| FCSA00141822.fits.gz | 2013-06-07 | DOMEFLAT | 11:26:05.681 | +19:54:09.65 | DOMEFLAT | 1.1 | 1.0 | 89.9361 | SCFCGRHD90 | 0.5 | -90.3 | SCFCSLLC20 | NONE | SCFCFLSO58 | NONE | 04:45:27.990 | 04:45:29.664 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| FCSA00141730.fits.gz | 2013-06-06 | BIAS | 23:03:44.623 | -01:20:01.29 | BIAS | 0.0 | 1.169 | 58.8068 | SCFCGRHD90 | 0.5 | -0.3 | SCFCSLLC20 | NONE | SCFCFLSO58 | NONE | 14:52:13.406 | 14:52:13.410 |
| FCSA00141818.fits.gz | 2013-06-07 | DOMEFLAT | 11:25:08.002 | +19:54:09.51 | DOMEFLAT | 1.1 | 1.0 | 89.9361 | SCFCGRHD90 | 0.5 | -90.3 | SCFCSLLC20 | NONE | SCFCFLSO58 | NONE | 04:44:30.393 | 04:44:31.697 |
| FCSA00142010.fits.gz | 2013-06-07 | Feige110 | 23:20:01.964 | -05:10:51.68 | OBJECT | 120.0 | 1.208 | 55.8712 | SCFCGRHD90 | 0.5 | -405.3 | SCFCSLLC20 | NONE | SCFCFLSO58 | NONE | 15:03:10.694 | 15:05:11.030 |
| FCSA00141798.fits.gz | 2013-06-07 | DOMEFLAT | 11:18:59.560 | +19:54:08.48 | DOMEFLAT | 3.0 | 1.0 | 89.9361 | SCFCGRHD90 | 0.5 | -90.3 | SCFCSLLC08 | NONE | SCFCFLSO58 | NONE | 04:38:23.184 | 04:38:26.574 |
| FCSA00141854.fits.gz | 2013-06-07 | BIAS | 11:41:06.075 | +19:54:11.27 | BIAS | 0.0 | 1.0 | 89.9361 | SCFCGRMR01 | 0.5 | -90.3 | SCFCSLLC20 | NONE | SCFCFLSO58 | NONE | 05:00:25.637 | 05:00:25.674 |
| FCSA00141586.fits.gz | 2013-06-06 | BIAS | 11:32:20.002 | +19:54:10.54 | BIAS | 0.0 | 1.0 | 89.9625 | SCFCGRHD90 | 0.5 | -90.3 | SCFCSLLC20 | NONE | SCFCFLSO58 | NONE | 04:55:43.880 | 04:55:43.883 |
| FCSA00141562.fits.gz | 2013-06-06 | DOMEFLAT | 11:22:06.397 | +19:54:09.10 | DOMEFLAT | 1.1 | 1.0 | 89.9625 | SCFCGRHD90 | 0.5 | -90.3 | SCFCSLLC20 | NONE | SCFCFLSO58 | NONE | 04:45:32.080 | 04:45:33.717 |
| FCSA00141572.fits.gz | 2013-06-06 | DOMEFLAT | 11:25:16.188 | +19:54:09.60 | DOMEFLAT | 3.0 | 1.0 | 89.9625 | SCFCGRHD90 | 0.5 | -90.3 | SCFCSLLC08 | NONE | SCFCFLSO58 | NONE | 04:48:41.382 | 04:48:44.710 |
| FCSA00141728.fits.gz | 2013-06-06 | BIAS | 23:19:58.381 | -05:11:14.89 | BIAS | 0.0 | 1.248 | 53.1982 | SCFCGRHD90 | 0.5 | -0.3 | SCFCSLLC20 | NONE | SCFCFLSO58 | NONE | 14:51:48.193 | 14:51:48.197 |
| FCSA00141800.fits.gz | 2013-06-07 | DOMEFLAT | 11:19:30.656 | +19:54:08.58 | DOMEFLAT | 3.0 | 1.0 | 89.9361 | SCFCGRHD90 | 0.5 | -90.3 | SCFCSLLC08 | NONE | SCFCFLSO58 | NONE | 04:38:54.202 | 04:38:57.612 |
biastab = summary_table[summary_table["type"] == "BIAS"]
flattab = summary_table[summary_table["type"] == "DOMEFLAT"]
arctab = summary_table[
(summary_table["type"] == "COMPARISON") & (summary_table["slit"] == "SCFCSLLC08")
]
scitab = summary_table[summary_table["object"] == "GRB 130606A"]
stdtab_f110 = summary_table[
(summary_table["object"] == "Feige110") & (summary_table["type"] == "OBJECT")
]
stdtab_f34 = summary_table[
(summary_table["object"] == "Feige34") & (summary_table["type"] == "OBJECT")
]
# for displaying in zscale
interval = ZScaleInterval()
# plotting the images
fig, axes = plt.subplots(2, 2, figsize=(15, 4))
titles = ["Bias", "Flat", "Arc", "Sci"]
bias = fits.getdata(RAWDIR / biastab["filename"][0])
flat = fits.getdata(RAWDIR / flattab["filename"][0])
arc = fits.getdata(RAWDIR / arctab["filename"][0])
sci = fits.getdata(RAWDIR / scitab["filename"][0])
for i, img in enumerate([bias, flat, arc, sci]):
vmin, vmax = interval.get_limits(img)
ax = axes[i // 2][i % 2]
ax.imshow(img.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.text(
0.05,
0.95,
titles[i],
transform=ax.transAxes,
ha="left",
va="top",
bbox=dict(facecolor="white", alpha=0.8),
)
ax.tick_params(axis="both", length=0.0, labelleft=False, labelbottom=False)
plt.tight_layout()
del bias, flat, arc, sci
Preprocessing of the Raw Images#
Preprocessing of the raw images is the first step in the data reduction process. The raw images are usually not ready for scientific analysis and need to be preprocessed before they can be used for scientific analysis. The preprocessing steps include the following:
Bias subtraction
Dark subtraction
Flat fielding
Cosmic ray removal
There are lots of software available for preprocessing the raw images. In this tutorial, we will use the ccdproc package for preprocessing the raw images. The ccdproc package is a Python package for reducing astronomical data taken with CCD cameras. It provides many of the necessary tools for processing of CCD images. For more information on the ccdproc package, please visit the documentation, and for the complete guide on the reduction process for astronomical images using Python, CCD Reduction and Photometry Guide.
In this notebook, we will preprocess the raw images using preproc.py script. This script is available in the tutorial directory. Dark subtraction will not be performed in this notebook since the amount of dark current in the images is negligible.
Making the Master Bias#
Before making the master bias, we need to check the bias images for any defects. We skip this step here, but it is important to check since the bias images are used to calibrate all the other images. If there are defects in the bias images, they will be propagated to all the other images.
The master bias is made by combining all the bias images. The master bias will be subtracted from all the images to remove the bias level. For combining the bias images, the sigma-clipping algorithm is used to reject outliers and then the average of the remaining images is calculated.
biaslist = [RAWDIR / f for f in biastab["filename"]]
_ = preproc.combine_bias(biaslist, outdir=OUTDIR, outname="MBias.fits")
del _
Show code cell output
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
WARNING: FITSFixedWarning: PC001001= 1.00000000 / Pixel Coordinate translation matrix
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC001002= 0.00000000 / Pixel Coordinate translation matrix
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC002001= 0.00000000 / Pixel Coordinate translation matrix
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC002002= 1.00000000 / Pixel Coordinate translation matrix
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / The equatorial coordinate system
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T04:59:10.645' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'unitfix' made the change 'Changed units:
'degree' -> 'deg'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'celfix' made the change 'Unmatched celestial axes'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T14:53:03.880' from MJD-END'. [astropy.wcs.wcs]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T04:56:58.652' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T04:59:35.780' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T15:06:35.545' from MJD-END'. [astropy.wcs.wcs]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T05:00:00.751' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T04:56:34.054' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T15:07:00.686' from MJD-END'. [astropy.wcs.wcs]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T15:07:53.456' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T04:56:09.042' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T14:53:29.035' from MJD-END'. [astropy.wcs.wcs]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T15:06:10.558' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T04:55:19.314' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T04:58:45.726' from MJD-END'. [astropy.wcs.wcs]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T14:52:38.925' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T15:07:26.756' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T14:52:13.481' from MJD-END'. [astropy.wcs.wcs]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T05:00:25.744' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T04:55:43.954' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T14:51:48.267' from MJD-END'. [astropy.wcs.wcs]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
mbias = OUTDIR / "MBias.fits"
mbias_img = fits.getdata(mbias)
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
vmin, vmax = interval.get_limits(mbias_img)
im = ax.imshow(mbias_img.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Master Bias")
divider = make_axes_locatable(ax)
cax = divider.append_axes("bottom", size="5%", pad=0.5)
plt.colorbar(im, cax=cax, orientation="horizontal")
fig.show()
Making Master Flat#
Not like the photometric observation, spectroscopic flat image is typically taken with something like halogen lamp or tungsten lamp, which is inevitably not uniform especially for the dispersion direction. Therefore, we need to make a master flat image flat, which involves modeling the flat spectrum and dividing the flat image by the model. To obtain a good master flat, we need to:
Examine the flat spectra: do they show any variablity throughout the observation? can we combine them without any further consideration? do they have any features like wiggles or bumps? do they have any discontinuity?
Model the flat spectra: we can use a polynomial (or smoothing convolution) to model the flat spectrum. The degree of the polynomial (or size of the smoothing kernel) should be determined by the complexity of the flat spectrum.
Normalize the flat spectra: divide the flat image by the model to obtain a normalized flat image.
Combining can be done either before or after the modeling. If the flat spectra are stable, we can combine them first and then model the combined flat spectrum. If the flat spectra are not stable, we can model each flat spectrum and then combine the normalized flat images. If the flat spectra show significant variability, we can also consider not combining them at all.
Because of this complexity, making a master flat can be a time-consuming part of the data reduction process. Moreover, the quality of the master flat directly affects the quality of the reduced science images, which may introduce systematic errors if not done properly. Therefore, the spectroscopic flat fielding is often omitted in the data reduction process, when the flat fielding conversely introduces more errors than it corrects.
Checking for Flat Image Shift#
The flat image may shift slightly depending on the telescope’s pointing direction. Since different images might require different flats, it can be risky to median combine all the flat images taken overnight. We should perform a quick and rough check to see if there is a shift in the flat images.
To avoid bias might arise from this shift, it is recommended to take flat images at the same pointing direction as the science images. However, if the shift is small, we can still use the median combined flat image for calibration.
Note
With the same reason, we should be careful when we reduce the arc frames. The location of the lines in the arc spectrum may shift slightly depending on the telescope’s pointing direction. Therefore, it is strongly recommended to take arc frames at the same pointing direction as the science frames.
We should check for the shift in the arc frames, but we will skip this process for simplicity.
However, unlike the flat fielding, the wavelength calibration is of critical importance in the spectroscopic data reduction process, since few percent of errors in the wavelength calibration can introduce systematic errors in the scientific analysis, such as the measurement of the redshift of galaxies. Therefore, we should always check for the shift in the arc frames before wavelength calibration.
In this notebook, we will skip this process for simplicity, but you should always check for flat image shift in your own data. We will assume that the flat images are already aligned, and we will proceed to combine them.
flattab_exp3 = flattab[flattab["exptime"] == 3.0]
flattab_exp1 = flattab[flattab["exptime"] == 1.1]
flatlist_exp3 = [RAWDIR / f for f in flattab_exp3["filename"]]
flatlist_exp1 = [RAWDIR / f for f in flattab_exp1["filename"]]
mbias = OUTDIR / "MBias.fits"
_ = preproc.make_master_flat(
flatlist_exp3,
outdir=OUTDIR,
outname="MFlat03.0.fits",
mbias=mbias,
filter_key="FILTER02",
verbose=False,
)
_ = preproc.make_master_flat(
flatlist_exp1,
outdir=OUTDIR,
outname="MFlat01.1.fits",
mbias=mbias,
filter_key="FILTER02",
verbose=False,
)
del _
Comparing with the Standard Star Profile#
Typical flat spectra shows many wiggles and bumps. Yet, we don’t know if these features are originated from the grating (as a representation of common instruments that can affect the both flat and science frames), which in turn will affect our science frame equally. If they are not due to the grating, such as filters in front of the flat-field lamp, the extreme difference in color temperature between the lamps and the object of interest, or wavelength-dependent variations in the reflectivity of the paint used for the flat-field screen (Massey & Hansen, 2010), they are not going to present on the science frame.
Let’s see if a supposedly smooth source (in this case, a standard star) has same “wiggles” and “bumps” on its spectrum. If this is the case, we should normalize the flat spectrum with a constant value (or a “low-order” functional model). If not, we should adopt “high-order” model or constant to normalize our flat field. For more information, please refer to Section 3.1., Massey & Hansen (2010).
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
flat03 = fits.getdata(OUTDIR / "MFlat03.0.fits")
vmin, vmax = interval.get_limits(flat03)
im = ax.imshow(flat03.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Master Flat (3.0s)")
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
flat01 = fits.getdata(OUTDIR / "MFlat01.1.fits")
vmin, vmax = interval.get_limits(flat01)
im = ax.imshow(flat03.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Master Flat (1.1s)")
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
f34 = fits.getdata(RAWDIR / stdtab_f34["filename"][0])
vmin, vmax = interval.get_limits(f34)
im = ax.imshow(f34.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Feige34 (Standard Star)")
Text(0.5, 1.0, 'Feige34 (Standard Star)')
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(211)
ax.plot(flat03[:, 700])
ax.text(
0.05,
0.95,
"Master Flat, 3.0s, Row 700",
transform=ax.transAxes,
ha="left",
va="top",
)
ax.grid()
ax = fig.add_subplot(212)
ax.plot(f34[:, 737])
ax.text(
0.05, 0.95, "Feige34, Row 700, Raw", transform=ax.transAxes, ha="left", va="top"
)
ax.grid()
Although we can see some tiny wiggles in the flat spectrum, this is not significant enough to adopt higher-order polynomials to model the flat spectrum, as opposed to the case if we have significant fringe patterns. If you closely examine the flat spectrum, you will see that the wiggles are not random, but rather located at the atmospheric absorption lines. This is because the flat image is taken with a lamp, and the lamp light is absorbed by the atmosphere at these absorption lines. Since the absorption lines are not random, we can safely ignore them when modeling the flat spectrum.
Here we model the flat spectrum with median filtering, which is a simple and effective way to remove noise while preserving the shape of the spectrum. We will use a median filter in the dispersion direction to model the flat spectrum, with a window size of 100 pixels. This is a good starting point, but you may need to adjust the window size depending on the complexity of the flat spectrum.
Note
Although I compared the horizontal profile of the standard star and the flat image for convenience, it might as well natural to compare the final extracted spectra of the standard star after applying both the high-order and low-order model of the flat field.
For simplicity, we will only use the combined flat image taken with the exposure time of 3 seconds. In practice, you may use all the flat images taken with different exposure times to make a master flat, or you can use the flat images taken with long and short exposure times to make a bad pixel map.
# median filtering the flat field image along the dispersion axis
flat03_filt = median_filter(flat03, size=50, axes=0)
# normalizing the flat field image
normflat = flat03 / flat03_filt
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(211)
ax.plot(flat03[:, 700], c="k", label="Raw", lw=1)
ax.plot(flat03_filt[:, 700], c="r", label="Filtered", lw=0.8)
ax.set_xlim(500, 4100)
ax.set_ylim(0, 2.2)
ax.legend()
ax.grid()
ax = fig.add_subplot(212)
ax.plot(normflat[:, 700], c="k", label="Normalized", lw=1)
ax.set_xlim(500, 4100)
ax.set_ylim(0.95, 1.05)
ax.legend()
ax.grid()
fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(311)
vmin, vmax = interval.get_limits(flat03)
im = ax.imshow(flat03.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Master Flat")
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="2%", pad=0.05)
plt.colorbar(im, cax=cax)
ax = fig.add_subplot(312)
im = ax.imshow(flat03_filt.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Filtered Master Flat")
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="2%", pad=0.05)
plt.colorbar(im, cax=cax)
ax = fig.add_subplot(313)
vmin, vmax = interval.get_limits(normflat)
im = ax.imshow(normflat.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Normalized Master Flat")
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="2%", pad=0.05)
plt.colorbar(im, cax=cax)
fig.show()
# writing the normalized flat field image to a file
hdr = fits.getheader(OUTDIR / "MFlat03.0.fits")
hdu = fits.PrimaryHDU(normflat, header=hdr)
hdu.writeto(OUTDIR / "nMFlat03.0.fits", overwrite=True)
Preprocessing the Raw Images#
Now we will preprocess the raw images. We will subtract the master bias and divide the master flat from the raw science frames. These preprocessed science frames will not be combined, since the spatial location of the object is different in each frame. We will combine the light from the object after extracting the 1D spectrum.
mbias = OUTDIR / "MBias.fits"
mflat = OUTDIR / "nMFlat03.0.fits"
# preprocessing the arc image
arc_list = [RAWDIR / f for f in arctab["filename"]]
print("Preprocessing arc images...")
arc = preproc.preproc(arc_list, mbias=mbias, mflat=mflat, outdir=OUTDIR)
# combine the preprocessed arc image
print("Combining the preprocessed arc images...")
parcs = [CCDData.read(OUTDIR / ("p" + f.stem)) for f in arc_list]
carc = combine(
parcs,
method="average",
sigma_clip=True,
sigma_clip_low_thresh=3,
sigma_clip_high_thresh=3,
)
carc.write(OUTDIR / "cArc.fits", overwrite=True)
print("The combined arc image is saved as cArc.fits")
# preprocessing the science image
print("Preprocessing science images...")
sci_list = [RAWDIR / f for f in scitab["filename"]]
_ = preproc.preproc(sci_list, mbias=mbias, mflat=mflat, outdir=OUTDIR, insert_ivar=True)
# preprocessing the standard star images
print("Preprocessing standard star images...")
print("Standard star: Feige110")
std_list_f110 = [RAWDIR / f for f in stdtab_f110["filename"]]
_ = preproc.preproc(
std_list_f110, mbias=mbias, mflat=mflat, outdir=OUTDIR, insert_ivar=True
)
print("Standard star: Feige34")
std_list_f34 = [RAWDIR / f for f in stdtab_f34["filename"]]
_ = preproc.preproc(
std_list_f34, mbias=mbias, mflat=mflat, outdir=OUTDIR, insert_ivar=True
)
Show code cell output
Preprocessing arc images...
Done: pFCSA00141840.fits
/Users/hbahk/class/TAOtutorials/tutorials/preproc.py:295: RuntimeWarning: invalid value encountered in divide
sci_data /= mflat_data # Flat fielding
Done: pFCSA00141882.fits
Done: pFCSA00141580.fits
Combining the preprocessed arc images...
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
WARNING: FITSFixedWarning: PC001001= 1.00000000 / Pixel Coordinate translation matrix
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC001002= 0.00000000 / Pixel Coordinate translation matrix
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC002001= 0.00000000 / Pixel Coordinate translation matrix
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: PC002002= 1.00000000 / Pixel Coordinate translation matrix
this form of the PCi_ja keyword is deprecated, use PCi_ja. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / The equatorial coordinate system
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T04:55:27.734' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'unitfix' made the change 'Changed units:
'degree' -> 'deg'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'celfix' made the change 'Unmatched celestial axes'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56450.000000 from DATE-OBS.
Set DATE-END to '2013-06-07T05:51:12.249' from MJD-END'. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56449.000000 from DATE-OBS.
Set DATE-END to '2013-06-06T04:52:48.140' from MJD-END'. [astropy.wcs.wcs]
INFO: An exception happened while extracting WCS information from the Header.
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3140 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.nddata.ccddata]
The combined arc image is saved as cArc.fits
Preprocessing science images...
Done: pFCSA00141920.fits
/Users/hbahk/class/TAOtutorials/tutorials/preproc.py:295: RuntimeWarning: invalid value encountered in divide
sci_data /= mflat_data # Flat fielding
Done: pFCSA00141930.fits
Done: pFCSA00141916.fits
Done: pFCSA00141928.fits
Done: pFCSA00141938.fits
Done: pFCSA00141940.fits
Done: pFCSA00141918.fits
Done: pFCSA00141936.fits
Done: pFCSA00141926.fits
Preprocessing standard star images...
Standard star: Feige110
Done: pFCSA00142004.fits
Done: pFCSA00142006.fits
Done: pFCSA00142010.fits
Standard star: Feige34
Done: pFCSA00141880.fits
Done: pFCSA00141894.fits
Done: pFCSA00141878.fits
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
vmin, vmax = interval.get_limits(carc.data)
im = ax.imshow(carc.data.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Combined Image")
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
sci = fits.getdata(OUTDIR / f"p{sci_list[0].stem}")
vmin, vmax = interval.get_limits(sci)
im = ax.imshow(sci.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Preprocessed Science Image")
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
std34 = fits.getdata(OUTDIR / f"p{std_list_f34[0].stem}")
vmin, vmax = interval.get_limits(std34)
im = ax.imshow(std34.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Preprocessed Feige34 Image")
fig.show()
Extracting 1D Spectra#
The 1D spectra are extracted from the 2D science frames by summing the flux along the spatial direction. The extraction process consists of the following steps:
Source Identification: Identify the source in the 2D science frame.
Aperture Tracing: Trace the spatial position of the source along the dispersion direction.
Aperture Summation: Sum the flux along the trace to extract the 1D spectrum.
Let’s extract the 1D spectra of the standard star first, because it is brighter and easier to extract.
1. Source Identification#
Let’s trim the images to the region of interest. Then, after the cosmic ray removal, we will identify the source in the 2D science frame. We will use the peak_local_max function from skimage.feature module to identify the source in a slice of the 2D science frame. The peak_local_max function finds the local maxima, which can be used to identify the source in the 2D science frame.
# Trim the images
TRIM = [0, 4220, 670, 810] # y_lower, y_upper, x_lower, x_upper
fname = OUTDIR / f"p{std_list_f34[0].stem}"
stdhdu = fits.open(fname)
stdhdr = stdhdu[0].header
stdtrim = stdhdu[0].data[TRIM[0] : TRIM[1], TRIM[2] : TRIM[3]]
exptime = stdhdr["EXPTIME"]
rdnoise = stdhdr["RDNOISE"]
fig = plt.figure(figsize=(12, 2))
ax = fig.add_subplot(111)
vmin, vmax = interval.get_limits(stdtrim)
im = ax.imshow(stdtrim.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Preprocessed Feige34 Image (Trimmed)")
fig.show()
The TRIM parameter is used to trim the images to the region of interest. The elements of the TRIM parameter is determined by examining the 2D science frame, to make the source centered in the image.
Cosmic ray removal often takes a long time. Hence, we will apply cosmic ray removal on the trimmed images to save unnecessary computation. To get reliable removal of cosmic rays, we need to mask bad pixels, but here we will skip this step for simplicity. For more information on cosmic ray removal, please refer to the CCD Data Reduction Guide.
# Cosmic ray removal
stdccd = CCDData(data=stdtrim, header=stdhdr, unit="adu")
gcorr = gain_correct(stdccd, gain=stdhdr["GAIN"] * u.electron / u.adu)
crrej = cosmicray_lacosmic(
gcorr, sigclip=7, readnoise=stdhdr["RDNOISE"] * u.electron, verbose=True
)
Starting 4 L.A.Cosmic iterations
Iteration 1:
325 cosmic pixels this iteration
Iteration 2:
34 cosmic pixels this iteration
Iteration 3:
21 cosmic pixels this iteration
Iteration 4:
19 cosmic pixels this iteration
Show code cell source
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(111)
vmin, vmax = interval.get_limits(crrej.data)
im = ax.imshow(crrej.data.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax)
ax.set_title("Cosmic Ray Removed Feige34 Image")
fig.show()
Looks good. Now we will identify the source position.
# Slice the image in the spatial direction
lcut = 1000 # left cut (dispersion direction)
rcut = 1050 # right cut (dispersion direction)
apall_probe = np.sum(crrej.data[lcut:rcut, :], axis=0)
x_probe = np.arange(apall_probe.size)
peak_pix = peak_local_max(
apall_probe, num_peaks=1, min_distance=50, threshold_abs=np.median(apall_probe)
)[0][0]
print("Peak pixel:", peak_pix)
Peak pixel: 68
Show code cell source
fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot(111)
ax.step(x_probe, apall_probe, where="mid")
ax.plot(
(peak_pix, peak_pix),
(
apall_probe[peak_pix] + 0.01 * apall_probe[peak_pix],
apall_probe[peak_pix] + 0.05 * apall_probe[peak_pix],
),
c="tab:red",
)
ax.annotate(
peak_pix,
(peak_pix, apall_probe[peak_pix] + 0.05 * apall_probe[peak_pix]),
textcoords="offset points",
xytext=(0, 10),
ha="left",
va="top",
color="k",
fontsize="small",
rotation=45,
)
ax.set_title("Aperture Sum")
ax.set_xlabel("Pixel Number (Spatial Direction)")
ax.set_ylabel("Counts [e-]")
fig.show()
Here we crop the image in the dispersion direction, from lcut to rcut. We sum the counts in the dispersion direction to get the 1D spatial profile of the source. We then use the peak_local_max function to identify the source position in the spatial profile. The peak_local_max function finds the local maxima in the 1D profile, which can be used to identify the source position. Then we can also estimate the full width at half maximum (FWHM) of the source, using the peak_widths function from the scipy.signal module.
From the above figure, we were able to get the pixel position of the peak of the source. However, the peak position is not always the center of the source, since it is an integer value. Therefore, we need some centering algorithm to get the center of the source to properly trace the spatial position of the source along the dispersion axis.
Here we use the centroid_com function from the photutils package to get the center of the source. The centroid_com function calculates the center of mass of the source, which is a good approximation of the center of the source. This function originally calculate the centroid of the source in the 2D image, but we will use them to calculate the center of the source in the 1D spatial profile. For more information on centroiding, please refer to the photutils documentation.
Before we estimate the center of the source, we need to estimate the background level. The background level is estimated by fitting a polynomial to the 1D spatial profile of the source. The background level is then subtracted from the 1D spatial profile to get the background-subtracted profile. The background-subtracted profile is then used to estimate the center of the source.
Let’s estimate the background level and subtract it from the 1D spatial profile. Here we use a Chebyshev polynomial of degree 3. The degree of the polynomial can be adjusted depending on the complexity of the background. To model the background, we should set the sky region where the source is not present. The sky region should be chosen carefully to avoid any contamination from the source.
# Select the sky region for the background estimation
sky_limit_peak = 10 # pixel limit from the peak pixel
mask_sky = np.ones_like(apall_probe, dtype=bool)
mask_sky[peak_pix - sky_limit_peak : peak_pix + sky_limit_peak] = False
x_probe_sky = x_probe[mask_sky]
apall_probe_sky = apall_probe[mask_sky]
fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot(111)
ax.step(x_probe, apall_probe, where="mid")
ax.fill_between(x_probe, 0, 1, where=mask_sky, alpha=0.1, transform=ax.get_xaxis_transform())
ax.set_xlabel("Pixel Number (Spatial Direction)")
ax.set_ylabel("Counts [e-]")
Text(0, 0.5, 'Counts [e-]')
# Fit the sky background with a Chebyshev polynomial
# define function for background fitting
def fit_background(x, apall, fitmask, sigma_sigclip, order_skymodel):
x_sky = x[fitmask]
apall_sky = apall[fitmask]
# sigma clipping
mask_sigclip = sigma_clip(apall_sky, sigma=sigma_sigclip, masked=True).mask
# fitting the sky background
coeffs_sky, fitfull_sky = chebfit(
x_sky[~mask_sigclip],
apall_sky[~mask_sigclip],
deg=order_skymodel,
full=True,
)
# sky background model
x = np.arange(apall.size)
skymodel = chebval(x, coeffs_sky)
# calculating the rms
residual_skymodel = fitfull_sky[0][0]
rms_skymodel = np.sqrt(residual_skymodel / len(x_sky[~mask_sigclip]))
return skymodel, rms_skymodel, mask_sigclip
# setting the parameters
sigma_sigclip = 3
order_skymodel = 3
# fitting the sky background
skymodel, rms_skymodel, mask_sigclip = fit_background(
x_probe,
apall_probe,
mask_sky,
sigma_sigclip=sigma_sigclip,
order_skymodel=order_skymodel,
)
# number of pixels
npix_sky = len(x_probe_sky)
npix_skymodel = np.count_nonzero(~mask_sigclip)
Show code cell source
fig, axes = plt.subplots(2, 1, figsize=(12, 8), gridspec_kw={"height_ratios": [4, 1]})
ax = axes[0]
ax.step(x_probe, apall_probe, where="mid")
ax.fill_between(
x_probe, 0, 1, where=mask_sky, alpha=0.1, transform=ax.get_xaxis_transform()
)
ax.plot(
x_probe,
skymodel,
c="tab:red",
lw=1,
label=f"Chebyshev Sky Model ({npix_skymodel:d}/{npix_sky:d})",
)
ax.plot(
x_probe_sky[mask_sigclip],
apall_probe_sky[mask_sigclip],
"x",
c="tab:orange",
label=f"Sigma-clipped ({sigma_sigclip:.1f}-sigma)",
)
resax = axes[1]
resax.plot([0, x_probe[-1]], [0, 0], c="k", ls="--", lw=0.8)
resax.step(
x_probe, apall_probe - skymodel, where="mid", label=f"RMS={rms_skymodel:.2f}"
)
resax.fill_between(
x_probe,
0,
1,
where=mask_sky,
color="tab:blue",
alpha=0.1,
transform=resax.get_xaxis_transform(),
)
resax.plot(
x_probe_sky[mask_sigclip],
apall_probe_sky[mask_sigclip] - skymodel[mask_sky][mask_sigclip],
"x",
c="tab:orange",
)
resax.set_xlabel("Pixel Number (Spatial Direction)")
resax.set_ylabel("Residuals [e-]")
# ylim_ax = ax.get_ylim()
ylim_resax = rms_skymodel * np.array([-1, 1]) * 10
resax.set_ylim(ylim_resax)
resax.tick_params(axis="x", direction="in", top="on")
resax.legend()
title_str = f"Skyfit: Chebyshev order {order_skymodel:d} ({sigma_sigclip:.1f}-sigma)"
ax.tick_params(axis="x", direction="in", labelbottom=False)
ax.set_xlabel("Pixel Number (Spatial Direction)")
ax.set_ylabel("Counts [e-]")
ax.set_title(title_str)
ax.legend()
for aa in axes:
aa.set_xlim(0, x_probe[-1])
fig.subplots_adjust(hspace=0)
From the sigma clipping, we can see some rejected points near the source. This indicates that the sky region we set is not sufficiently far from the source. We should set the sky region further away from the source to avoid any contamination from the source. Sigma clipping will also reject any cosmic rays or bad pixels in the sky region, which will affect the background estimation.
from photutils.centroids import centroid_com
apall_probe_skysub = apall_probe - skymodel
center_probe_com = centroid_com(apall_probe_skysub)[0]
print("Center without sky subtraction:", centroid_com(apall_probe)[0])
print("Center after sky subtraction:", center_probe_com)
pwresult = peak_widths(apall_probe_skysub, [peak_pix], rel_height=0.5)
fwhm_probe = pwresult[0][0]
left_ips, right_ips = pwresult[2][0], pwresult[3][0]
fwhm_hight = pwresult[1][0]
print("FWHM:", fwhm_probe)
print("Left and right interpolated positions:", left_ips, right_ips)
Center without sky subtraction: 67.87118517140294
Center after sky subtraction: 67.77667632982664
FWHM: 2.9532888855505917
Left and right interpolated positions: 66.33374812839358 69.28703701394417
Show code cell source
import plotly.io as pio
import plotly.express as px
import plotly.offline as py
import plotly.graph_objects as go
pio.renderers.default = "notebook"
fig = px.line(
x=np.arange(len(apall_probe_skysub)),
y=apall_probe,
title="Sky Subtracted Sum",
labels={"x": "Pixel Number (Spatial Direction)", "y": "Counts [e-]"},
line_shape="hvh"
)
fig.add_vline(
x=peak_pix,
line_dash="dash",
line_color="red",
annotation_text="Peak",
annotation_position="top right",
line_width=0.8,
)
fig.add_vline(
x=center_probe_com,
line_dash="dash",
line_color="green",
annotation_text="Center",
annotation_position="top left",
line_width=0.8,
)
fig.add_trace(go.Scatter(
x=[left_ips, right_ips],
y=[fwhm_hight, fwhm_hight],
mode="lines+markers",
line=dict(color="red", width=1.5),
name=f"FWHM={fwhm_probe:.2f}",
))
fig
Zoom in to see the subtle difference between the peak position and the center position.
2. Aperture Tracing#
We have identified the source position in the spatial direction by collapsing a slice of the 2D science frame along the dispersion direction. Next, we will trace the spatial position of the source along the dispersion direction. This involves changing the location of the slices in the dispersion direction and identifying the source position in each slice.
# Define the function for determining the center of the source for given slice range
def get_skymask(apall, sky_limit_peak):
peak_pix = peak_local_max(
apall,
num_peaks=1,
min_distance=10,
threshold_abs=np.mean(apall),
)[0][0]
mask_sky = np.ones_like(apall, dtype=bool)
mask_sky[peak_pix - sky_limit_peak : peak_pix + sky_limit_peak] = False
return mask_sky, peak_pix
def find_center(data, lcut, rcut, sky_limit_peak, sigma_sigclip, order_skymodel):
apall = np.sum(data[lcut:rcut, :], axis=0)
x = np.arange(apall.size)
# mask for the sky region and initial peak pixel
mask_sky, peak_pix = get_skymask(apall, sky_limit_peak)
# fitting the sky background
skymodel, rms_skymodel, mask_sigclip = fit_background(
x,
apall,
mask_sky,
sigma_sigclip=sigma_sigclip,
order_skymodel=order_skymodel,
)
apall_skysub = apall - skymodel
# center of the source
center_com = centroid_com(apall_skysub)[0]
# fwhm of the source
pwresult = peak_widths(apall_skysub, [round(peak_pix)], rel_height=0.5)
fwhm = pwresult[0][0]
return center_com, fwhm
# Find the center of the source for all the slices
slice_width = 50
sky_limit_peak = 20 # increased the distance from the source
midpoints = []
aptrace = []
aptrace_fwhm = []
lcut, rcut = 0, slice_width
while rcut < len(crrej.data):
rcut = lcut + slice_width
midpoints.append((lcut + rcut) / 2)
center, fwhm = find_center(
crrej.data, lcut, rcut, sky_limit_peak, sigma_sigclip, order_skymodel
)
aptrace.append(center)
aptrace_fwhm.append(fwhm)
lcut += slice_width
rcut += slice_width
aptrace = np.array(aptrace)
aptrace_fwhm = np.array(aptrace_fwhm)
aptrace_sigma = aptrace_fwhm * gaussian_fwhm_to_sigma
midpoints = np.array(midpoints)
# Let's fit the aperture trace
order_aptrace = 9
sigma_clip_aptrace = 3
# fitting the trace line
coeffs_aptrace_init = chebfit(midpoints, aptrace, deg=order_aptrace)
# sigma clipping
mask_aptrace = sigma_clip(
aptrace - chebval(midpoints, coeffs_aptrace_init),
sigma=sigma_clip_aptrace,
masked=True,
).mask
# fitting the trace line again
coeffs_aptrace = chebfit(
midpoints[~mask_aptrace], aptrace[~mask_aptrace], deg=order_aptrace
)
residuals_aptrace = aptrace - chebval(midpoints, coeffs_aptrace)
rejected = ~np.in1d(midpoints, midpoints[~mask_aptrace])
rms_aptrace = np.sqrt(np.mean(residuals_aptrace[~mask_aptrace] ** 2))
# aperture trace model
ximg = np.arange(crrej.data.shape[0])
aptrace_model = chebval(ximg, coeffs_aptrace)
Show code cell source
# plot the aperture trace
fig, axes = plt.subplots(2, 1, figsize=(12, 6), gridspec_kw={"height_ratios": [3, 1]})
ax = axes[0]
ax.imshow(
crrej.data.T, cmap="gray", origin="lower", vmin=vmin, vmax=vmax, aspect="auto"
)
ax.scatter(midpoints, aptrace, marker="+", c="dodgerblue", s=100, lw=2, label="Data")
ax.scatter(
midpoints[rejected],
aptrace[rejected],
marker="x",
c="salmon",
s=100,
lw=1,
label="Rejected",
)
ax.plot(ximg, aptrace_model, c="tab:red", lw=1, label="Chebyshev Trace Model")
ax.plot(midpoints, aptrace + 5*aptrace_sigma, c="r", lw=0.8, ls="--", label="5-sigma")
ax.plot(midpoints, aptrace - 5*aptrace_sigma, c="r", lw=0.8, ls="--")
ax.plot(midpoints, aptrace + 10*aptrace_sigma, c="r", lw=0.8, ls=":", label="10-sigma")
ax.plot(midpoints, aptrace - 10*aptrace_sigma, c="r", lw=0.8, ls=":")
ax.set_title("Aperture Trace")
ax.set_ylabel("Pixel Number (Spatial Direction)")
ax.set_xticks([])
ax.legend()
resax = axes[1]
resax.plot([0, ximg[-1]], [0, 0], c="k", ls="--", lw=0.8)
resax.scatter(midpoints, residuals_aptrace, marker="+", c="dodgerblue", s=100, lw=2)
resax.scatter(
midpoints[rejected],
residuals_aptrace[rejected],
marker="x",
c="salmon",
s=100,
lw=1,
)
resax.set_xlabel("Pixel Number (Dispersion Direction)")
resax.set_ylabel("Residuals [pix]")
resax.tick_params(axis="x", direction="in", top="on")
resax_ylim = np.max(np.abs(residuals_aptrace)) * np.array([-1, 1])
resax.set_ylim(resax_ylim)
resax.set_xlim(0, ximg[-1])
resax.text(
0.99,
0.95,
f"RMS={rms_aptrace:.2f}",
transform=resax.transAxes,
ha="right",
va="top",
)
fig.subplots_adjust(hspace=0)
3. Aperture Summation#
Now we have traced the spatial position of the source along the dispersion direction. We will sum the flux along the trace to extract the 1D spectrum. The aperture size should be chosen carefully to avoid any contamination from the sky or other sources. The aperture size should be large enough to capture sufficient amount of the flux from the source but small enough to avoid any contamination from the sky or other sources.
See also
The type of aperture used for the extraction can affect the signal-to-noise ratio of the extracted spectrum. One can use a simple aperture, which sums the flux within a fixed aperture size, or an optimal extraction (Horne 1986), which weights the flux based on the spatial profile of the source. The optimal extraction can improve the signal-to-noise ratio of the extracted spectrum, especially when the source is extended or the background is not uniform. In this tutorial, we will use a simple aperture for the extraction, but you can try the optimal extraction for better results.
Let’s make a simple aperture summation to extract the 1D spectrum!
# Define the function for aperture summation
def aperture_sum(
cut, center, apsum_sigma, ap_sigma, sky_limit_peak, sigma_sigclip, order_skymodel
):
x = np.arange(len(cut))
mask_sky, peak_pix = get_skymask(cut, sky_limit_peak)
# aperture size = trace center +/- apsum_sigma * ap_sigma
appos_lower = int(center - apsum_sigma * ap_sigma)
appos_upper = int(center + apsum_sigma * ap_sigma)
# fitting the sky background
skymodel, rms_skymodel, mask_sigclip = fit_background(
x,
cut,
mask_sky,
sigma_sigclip=sigma_sigclip,
order_skymodel=order_skymodel,
)
# subtracting the sky background
cut_skysub = cut - skymodel
source = cut_skysub[appos_lower:appos_upper]
# aperture sum
ap_sum = np.sum(source)
# calculating the uncertainty
sig = rdnoise**2 + source + skymodel[appos_lower:appos_upper]
ap_sig = np.sqrt(np.sum(sig))
return ap_sum, ap_sig
# Parameters for the aperture summation
apsum_sigma = 10
ap_fwhm = np.median(aptrace_fwhm[~mask_aptrace]) # median FWHM in pixels
ap_sigma = ap_fwhm * gaussian_fwhm_to_sigma # sigma under the Gaussian assumption
# Extract the spectrum along the dispersion direction
ap_sum = []
ap_sig = []
for i in range(crrej.data.shape[0]):
cut_i = crrej.data[i, :]
center_i = aptrace_model[i]
ap_sum_i, ap_sig_i = aperture_sum(
cut_i,
center_i,
apsum_sigma,
ap_sigma,
sky_limit_peak,
sigma_sigclip,
order_skymodel,
)
ap_sum.append(ap_sum_i)
ap_sig.append(ap_sig_i)
ap_sum = np.array(ap_sum) / exptime
ap_sig = np.array(ap_sig) / exptime
Show code cell source
# Plot the extracted spectrum
fig = go.Figure()
fig.add_trace(go.Scatter(
x=ximg,
y=ap_sum,
error_y=dict(
type="data",
array=ap_sig,
visible=True,
color="dodgerblue",
thickness=1,
width=1,
),
mode="lines+markers",
marker=dict(color="black", size=3),
line=dict(color="black", width=0.8),
))
fig.update_layout(
title="Extracted Spectrum of Feige34",
xaxis_title="Pixel Number (Dispersion Direction)",
yaxis_title="Counts per Exposure Time[e-/s]",
showlegend=False
)
fig
# Let's make a class for the extraction of the spectrum
class InstrumentalSpectrum:
def __init__(
self,
fpath,
sky_limit_peak,
sigma_sigclip,
order_skymodel,
order_aptrace,
sigma_clip_aptrace,
slice_width_aptrace,
apsum_sigma,
trim=[0, 4220, 670, 810],
):
self.fname = fpath
self.TRIM = trim
self.SKY_LIMIT_PEAK = sky_limit_peak
self.SIGMA_SIGCLIP = sigma_sigclip
self.ORDER_SKYMODEL = order_skymodel
self.ORDER_APTRACE = order_aptrace
self.SIGMA_CLIP_APTRACE = sigma_clip_aptrace
self.SLICE_WIDTH_APTRACE = slice_width_aptrace
self.APSUM_SIGMA = apsum_sigma
self.hdu = fits.open(fpath)
self.data = self.hdu[0].data[
self.TRIM[0] : self.TRIM[1], self.TRIM[2] : self.TRIM[3]
]
self.hdr = self.hdu[0].header
self.EXPTIME = self.hdr["EXPTIME"]
self.RDNOISE = self.hdr["RDNOISE"]
self.ccd = CCDData(data=self.data, header=self.hdr, unit="adu")
self.gcorr = gain_correct(self.ccd, gain=self.hdr["GAIN"] * u.electron / u.adu)
self.crrej = cosmicray_lacosmic(
self.gcorr,
sigclip=7,
readnoise=self.hdr["RDNOISE"] * u.electron,
verbose=False,
)
self.set_aptrace()
self.set_apsum()
def set_aptrace(self):
slice_width = self.SLICE_WIDTH_APTRACE
lcut, rcut = 0, slice_width
midpoints = []
aptrace = []
aptrace_fwhm = []
while rcut < len(self.crrej.data):
rcut = lcut + slice_width
midpoints.append((lcut + rcut) / 2)
center, fwhm = find_center(
self.crrej.data,
lcut,
rcut,
self.SKY_LIMIT_PEAK,
self.SIGMA_SIGCLIP,
self.ORDER_SKYMODEL,
)
aptrace.append(center)
aptrace_fwhm.append(fwhm)
lcut += slice_width
rcut += slice_width
self.aptrace = np.array(aptrace)
self.aptrace_fwhm = np.array(aptrace_fwhm)
midpoints = np.array(midpoints)
# fitting the trace line
coeffs_aptrace_init = chebfit(midpoints, self.aptrace, deg=self.ORDER_APTRACE)
# sigma clipping
mask_aptrace = sigma_clip(
self.aptrace - chebval(midpoints, coeffs_aptrace_init),
sigma=self.SIGMA_CLIP_APTRACE,
masked=True,
).mask
self.mask_aptrace = mask_aptrace
# fitting the trace line again
self.coeffs_aptrace = chebfit(
midpoints[~mask_aptrace],
self.aptrace[~mask_aptrace],
deg=self.ORDER_APTRACE,
)
self.residuals_aptrace = self.aptrace - chebval(midpoints, self.coeffs_aptrace)
self.rejected_aptrace = ~np.in1d(midpoints, midpoints[~mask_aptrace])
self.rms_aptrace = np.sqrt(np.mean(self.residuals_aptrace[~mask_aptrace] ** 2))
# aperture trace model
ximg = np.arange(self.crrej.data.shape[0])
self.aptrace_model = chebval(ximg, self.coeffs_aptrace)
def set_apsum(self):
ap_sum = []
ap_sig = []
self.ap_fwhm = np.median(self.aptrace_fwhm[~self.mask_aptrace])
self.ap_sigma = self.ap_fwhm * gaussian_fwhm_to_sigma
for i in range(self.crrej.data.shape[0]):
cut_i = self.crrej.data[i, :]
center_i = self.aptrace_model[i]
ap_sum_i, ap_sig_i = aperture_sum(
cut_i,
center_i,
self.APSUM_SIGMA,
self.ap_sigma,
self.SKY_LIMIT_PEAK,
self.SIGMA_SIGCLIP,
self.ORDER_SKYMODEL,
)
ap_sum.append(ap_sum_i)
ap_sig.append(ap_sig_i)
self.ap_sum = np.array(ap_sum) / self.EXPTIME
self.ap_sig = np.array(ap_sig) / self.EXPTIME
std_insspecs = []
for std in std_list_f34:
fname = OUTDIR / f"p{std.stem}"
ins = InstrumentalSpectrum(
fname,
sky_limit_peak=20,
sigma_sigclip=3,
order_skymodel=3,
order_aptrace=9,
sigma_clip_aptrace=3,
slice_width_aptrace=50,
apsum_sigma=10,
)
std_insspecs.append(ins)
Show code cell source
fig = go.Figure()
f, axes = plt.subplots(3, 1, figsize=(12, 3))
i=0
for StdInsSpec in std_insspecs:
ax = axes[i]
fig.add_trace(go.Scatter(
x=ximg,
y=StdInsSpec.ap_sum,
error_y=dict(
type="data",
array=StdInsSpec.ap_sig,
visible=True,
thickness=1,
width=1,
),
mode="lines+markers",
marker=dict(size=3),
line=dict(width=0.8),
name=StdInsSpec.fname.stem,
))
i += 1
ax.imshow(
StdInsSpec.crrej.data.T,
cmap="gray",
origin="lower",
vmin=vmin,
vmax=vmax,
)
ax.set_title(StdInsSpec.fname.stem)
f.suptitle("2D Images of the Standard Star Frames")
fig.update_layout(
title="Extracted Spectrum of Feige34",
xaxis_title="Pixel Number (Dispersion Direction)",
yaxis_title="Counts per Exposure Time[e-/s]",
showlegend=True,
legend=dict(
yanchor="top",
y=0.99,
xanchor="left",
x=0.01
))
fig
stdtab_f34
| filename | date-obs | object | ra | dec | type | exptime | airmass | altitude | disperser | slit-width | slit-pa | slit | filter1 | filter2 | filter3 | ut-str | ut-end |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| str50 | str50 | str50 | str50 | str50 | str50 | float64 | float64 | float64 | str50 | float64 | float64 | str50 | str50 | str50 | str50 | str50 | str50 |
| FCSA00141880.fits.gz | 2013-06-07 | Feige34 | 10:39:30.972 | +43:06:49.26 | OBJECT | 60.0 | 1.188 | 57.2701 | SCFCGRHD90 | 0.5 | -240.3 | SCFCSLLC08 | NONE | SCFCFLSO58 | NONE | 05:49:08.496 | 05:50:08.794 |
| FCSA00141894.fits.gz | 2013-06-07 | Feige34 | 10:39:30.923 | +43:06:49.54 | OBJECT | 60.0 | 1.207 | 55.9261 | SCFCGRHD90 | 0.5 | -240.3 | SCFCSLLC20 | NONE | SCFCFLSO58 | NONE | 05:58:05.720 | 05:59:05.978 |
| FCSA00141878.fits.gz | 2013-06-07 | Feige34 | 10:39:30.978 | +43:06:49.28 | OBJECT | 60.0 | 1.185 | 57.5049 | SCFCGRHD90 | 0.5 | -240.3 | SCFCSLLC08 | NONE | SCFCFLSO58 | NONE | 05:47:32.380 | 05:48:32.720 |
As you can see from the figure above, the extracted spectrum of FCSA00142006 shows different level of continuum. It can be due to various reasons, such as slight shift of the source position with respect to the slit aperture, or the difference in the seeing condition, out-of-focus, etc.
To correct for this, we can normalize the extracted spectrum by dividing it by the continuum level at some range. Here we will normalize the extracted spectrum by dividing it by the median value of the continuum level at the range of 1000-2000 pixels.
Note
If we manipulate the continuum level of standard star like this, how should we calibrate the flux of the science frame and trust them? For this reason, we need a process called aperture correction to correct the flux. We will skip this process for simplicity, but you should always consider this process in your own data reduction.